traits
above- and belowground
acquisitive vs. conservative resource-use strategies
collaborative vs. individual resource-use strategies
fine root proportion
thickness of roots
Goal: Examine how traits shift in community compositional (ordinational) space, placing each trait on an axis.
Author: Caryn D. Iwanaga
Updated: 01/24/2025
library(tidyverse)
library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)
library(readxl)
cover.dat.labels <- read.csv("https://www.dropbox.com/scl/fi/up8nnzkcpchsr45f8cm92/Compost_Cover_LongClean.csv?rlkey=z2tvaj8t6khadef7ydz782zka&st=qwef9ys0&dl=1") %>%
mutate(
nut_trt = factor(ifelse(nut_trt =="C", "+compost", ifelse(nut_trt =="F", "+N fertilizer", ifelse(nut_trt =="N", "control", NA))), levels = c("control", "+N fertilizer", "+compost")),
ppt_trt = ifelse(ppt_trt == "D", "dry", ifelse(ppt_trt == "XC", "control", ifelse(ppt_trt == "W", "wet", NA))))
site.df <-
cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
mutate(
grazing_hist = ifelse(block == 1 | block == 2, "high", ifelse(block == 3 | block == 4, "low", NA))
) %>% # separate by block (grazing intensities)
group_by(nut_trt, ppt_trt, yr, code4, grazing_hist) %>% # make each code unique -- one value for species abundance for each trt
summarise(
# pct_cover,
abundance = mean(pct_cover)
# sd = sd(pct_cover)
)
`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.
matrix0 <- spread(
site.df,
code4,
abundance,
fill = 0)
# make rownames ----
rownames(matrix0) <- paste(matrix0$nut_trt, matrix0$ppt_trt, matrix0$yr, matrix0$grazing_hist, sep = "_")
Warning: Setting row names on a tibble is deprecated.
# delete columns ----
matrix <- matrix0[, c(5:79)]
# relativize data with Bray-Curtis dissimilarity
matrix.bray <- decostand(matrix, "total")
rownames(matrix.bray) <- rownames(matrix0)
site.block.sd.df <-
cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
summarise(
pct_cover,
abundance = mean(pct_cover),
sd = sd(pct_cover)
)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4', 'block'. You can override using the `.groups` argument.
Rules of thumb:
< 0.05 is excellent
0.05-0.1 is good
0.1-0.2 is fair
0.2-0.3 is cause for concern…
> 0.3 is poor, random
Stress: fair
0.1966553
Grouping by year and grazing history had most significant clustering shown.
nmds.comp <- metaMDS(matrix.bray, k=2, trymax = 25)
Run 0 stress 0.2020736
Run 1 stress 0.1970655
... New best solution
... Procrustes: rmse 0.06573517 max resid 0.2315863
Run 2 stress 0.2027495
Run 3 stress 0.1992376
Run 4 stress 0.2045952
Run 5 stress 0.2026252
Run 6 stress 0.2027495
Run 7 stress 0.2093589
Run 8 stress 0.2020736
Run 9 stress 0.1970625
... New best solution
... Procrustes: rmse 0.006290782 max resid 0.03755233
Run 10 stress 0.2075871
Run 11 stress 0.2133994
Run 12 stress 0.213313
Run 13 stress 0.2037336
Run 14 stress 0.2027496
Run 15 stress 0.2026295
Run 16 stress 0.2043011
Run 17 stress 0.20163
Run 18 stress 0.2059418
Run 19 stress 0.1973472
... Procrustes: rmse 0.02205658 max resid 0.1152837
Run 20 stress 0.2101999
Run 21 stress 0.1970631
... Procrustes: rmse 0.006111069 max resid 0.03747752
Run 22 stress 0.2171655
Run 23 stress 0.2017185
Run 24 stress 0.1974235
... Procrustes: rmse 0.02005767 max resid 0.1153592
Run 25 stress 0.2033605
*** Best solution was not repeated -- monoMDS stopping criteria:
23: stress ratio > sratmax
2: scale factor of the gradient < sfgrmin
nmds.comp # pretty not great stress
Call:
metaMDS(comm = matrix.bray, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray
Distance: bray
Dimensions: 2
Stress: 0.1970625
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 9 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray’
stressplot(nmds.comp)
plot(nmds.comp, type="t")
nmds1 <- as.data.frame(scores(nmds.comp, choices=c(1), display=c("sites"))) %>%
mutate(matrix0[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.comp, choices=c(2), display=c("sites")))
nmds_dat <- cbind(nmds1, nmds2) %>%
Error: Incomplete expression: select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year")
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Grazing History")
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(grazing_hist), fill = as.factor(yr))) +
geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "Grouped by Grazing History & Year") +
scale_color_manual(
values = c(
"low" = "cyan",
"high" = "brown",
"2019" = "salmon",
"2020" = "limegreen",
"2021" = "cornflowerblue")
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(nut_trt)),
geom = "polygon",
fill = NA, # Makes the ellipse transparent inside
size = 1) +
labs(title = "Grouped by Nutrient Treatment & Year") +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
Stress: cause for concern
0.2329306
nmds.block.comp <- metaMDS(matrix.block.bray, k=2, trymax = 25)
Run 0 stress 0.2329306
Run 1 stress 0.243463
Run 2 stress 0.2482815
Run 3 stress 0.2498999
Run 4 stress 0.2329786
... Procrustes: rmse 0.002230422 max resid 0.01280587
Run 5 stress 0.2333533
... Procrustes: rmse 0.01157363 max resid 0.1126004
Run 6 stress 0.2421491
Run 7 stress 0.2443351
Run 8 stress 0.2333565
... Procrustes: rmse 0.01146138 max resid 0.1123691
Run 9 stress 0.246737
Run 10 stress 0.2561412
Run 11 stress 0.2432135
Run 12 stress 0.2557287
Run 13 stress 0.2438602
Run 14 stress 0.2496649
Run 15 stress 0.2406222
Run 16 stress 0.2594277
Run 17 stress 0.2535404
Run 18 stress 0.2329445
... Procrustes: rmse 0.001578641 max resid 0.01195696
Run 19 stress 0.2329307
... Procrustes: rmse 4.0144e-05 max resid 0.000261745
... Similar to previous best
Run 20 stress 0.236599
*** Best solution repeated 1 times
nmds.block.comp
Call:
metaMDS(comm = matrix.block.bray, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.block.bray
Distance: bray
Dimensions: 2
Stress: 0.2329306
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.block.bray’
stressplot(nmds.block.comp)
plot(nmds.block.comp, type="t")
nmds1 <- as.data.frame(scores(nmds.block.comp, choices=c(1), display=c("sites"))) %>%
mutate(matrix0.block[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.block.comp, choices=c(2), display=c("sites")))
nmds_dat.block <- cbind(nmds1, nmds2) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment - Blocks") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year - Blocks")
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment - Blocks") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = as.factor(block), color = as.factor(block))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Block")
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = as.factor(block), fill = as.factor(yr))) +
geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(block)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "Grouped by Block & Year - Blocks") +
scale_color_manual(
values = c(
# "low" = "cyan",
# "high" = "brown",
"2019" = "salmon",
"2020" = "limegreen",
"2021" = "cornflowerblue")
)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(nut_trt)),
geom = "polygon",
fill = NA, # Makes the ellipse transparent inside
size = 1) +
labs(title = "Grouped by Nutrient Treatment & Year - Blocks") +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
Stress: fair
0.1869068
# relativize data with Bray-Curtis dissimilarity
matrix.low <- matrix0 %>%
filter(grazing_hist == "low")
row.low <- paste(matrix.low$nut_trt, matrix.low$ppt_trt, matrix.low$yr, matrix.low$grazing_hist, sep = "_")
matrix.low <- matrix.low[, c(5:79)]
rownames(matrix.low) <- row.low
Warning: Setting row names on a tibble is deprecated.
matrix.bray.low <- decostand(matrix.low, "total")
nmds.comp.low <- metaMDS(matrix.bray.low, k=2, trymax = 25)
Run 0 stress 0.2000574
Run 1 stress 0.1869068
... New best solution
... Procrustes: rmse 0.1206791 max resid 0.4645309
Run 2 stress 0.1941799
Run 3 stress 0.1876869
Run 4 stress 0.1930573
Run 5 stress 0.1869068
... Procrustes: rmse 8.072635e-05 max resid 0.0002281104
... Similar to previous best
Run 6 stress 0.1908791
Run 7 stress 0.1942011
Run 8 stress 0.2045135
Run 9 stress 0.2142868
Run 10 stress 0.2264618
Run 11 stress 0.1876848
Run 12 stress 0.201295
Run 13 stress 0.2013001
Run 14 stress 0.2298215
Run 15 stress 0.2303754
Run 16 stress 0.1941799
Run 17 stress 0.2066946
Run 18 stress 0.2237619
Run 19 stress 0.1941799
Run 20 stress 0.2122824
*** Best solution repeated 1 times
nmds.comp.low # pretty not great stress
Call:
metaMDS(comm = matrix.bray.low, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.low
Distance: bray
Dimensions: 2
Stress: 0.1869068
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 1 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.low’
stressplot(nmds.comp.low)
plot(nmds.comp.low, type="t")
nmds1 <- as.data.frame(scores(nmds.comp.low, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.low), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.low), "_"), `[`, 4)
)
nmds2 <- as.data.frame(scores(nmds.comp.low, choices=c(2), display=c("sites")))
nmds_low_dat <- cbind(nmds1, nmds2) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, color = as.factor(yr), shape = nut_trt)) +
geom_point(size = 3, stroke = 1) +
labs(title = "2019-2021 (LOW) - Grouping by Year & Nutrient Treatment")
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment")
# add year ellipses ----
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+
custom_colors
# relativize data with Bray-Curtis dissimilarity
matrix.2019 <- matrix0 %>%
filter(yr == 2019)
row.2019 <- paste(matrix.2019$nut_trt, matrix.2019$ppt_trt, matrix.2019$yr, matrix.2019$grazing_hist, sep = "_")
matrix.2019 <- matrix.2019[, c(5:79)]
rownames(matrix.2019) <- row.2019
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2019 <- decostand(matrix.2019, "total")
Stress: Excellent
0.0560178
nmds.2019 <- metaMDS(matrix.bray.2019, k=2, trymax = 25)
Run 0 stress 0.0560178
Run 1 stress 0.0560178
... Procrustes: rmse 2.857261e-06 max resid 9.034415e-06
... Similar to previous best
Run 2 stress 0.05905532
Run 3 stress 0.07335646
Run 4 stress 0.0560178
... Procrustes: rmse 5.903723e-06 max resid 1.375655e-05
... Similar to previous best
Run 5 stress 0.05960972
Run 6 stress 0.07030125
Run 7 stress 0.06384516
Run 8 stress 0.05877078
Run 9 stress 0.05960972
Run 10 stress 0.06223098
Run 11 stress 0.05877078
Run 12 stress 0.05654879
Run 13 stress 0.05905533
Run 14 stress 0.06677855
Run 15 stress 0.05877078
Run 16 stress 0.05960972
Run 17 stress 0.059968
Run 18 stress 0.05957995
Run 19 stress 0.05639216
... Procrustes: rmse 0.008347172 max resid 0.02870258
Run 20 stress 0.05639215
... Procrustes: rmse 0.008300352 max resid 0.02861537
*** Best solution repeated 2 times
nmds.2019 # pretty not great stress
Call:
metaMDS(comm = matrix.bray.2019, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2019
Distance: bray
Dimensions: 2
Stress: 0.0560178
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2019’
stressplot(nmds.2019)
plot(nmds.2019, type="t")
nmds1.2019 <- as.data.frame(scores(nmds.2019, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 4)
)
nmds2.2019 <- as.data.frame(scores(nmds.2019, choices=c(2), display=c("sites")))
nmds_2019_dat <- cbind(nmds1.2019, nmds2.2019) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2019_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 3, stroke = 1) +
labs(title = "2019 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2019_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
labs(title = "2019 - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_2019_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2019 - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2020 <- matrix0 %>%
filter(yr == 2020)
row.2020 <- paste(matrix.2020$nut_trt, matrix.2020$ppt_trt, matrix.2020$yr, matrix.2020$grazing_hist, sep = "_")
matrix.2020 <- matrix.2020[, c(5:79)]
rownames(matrix.2020) <- row.2020
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2020 <- decostand(matrix.2020, "total")
Stress: fair
0.1522495
nmds.2020 <- metaMDS(matrix.bray.2020, k=2, trymax = 25)
Run 0 stress 0.1524241
Run 1 stress 0.1748537
Run 2 stress 0.1584809
Run 3 stress 0.1748534
Run 4 stress 0.1692911
Run 5 stress 0.153917
Run 6 stress 0.1744239
Run 7 stress 0.1524241
... New best solution
... Procrustes: rmse 9.73555e-06 max resid 2.454011e-05
... Similar to previous best
Run 8 stress 0.1625903
Run 9 stress 0.1761789
Run 10 stress 0.1830193
Run 11 stress 0.1743562
Run 12 stress 0.1584809
Run 13 stress 0.1592892
Run 14 stress 0.1866315
Run 15 stress 0.1522045
... New best solution
... Procrustes: rmse 0.02481256 max resid 0.06590007
Run 16 stress 0.1524241
... Procrustes: rmse 0.02481195 max resid 0.06599272
Run 17 stress 0.1524241
... Procrustes: rmse 0.02480968 max resid 0.06599522
Run 18 stress 0.1633441
Run 19 stress 0.1594877
Run 20 stress 0.1522495
... Procrustes: rmse 0.01182949 max resid 0.0373637
Run 21 stress 0.1754284
Run 22 stress 0.1595535
Run 23 stress 0.1592892
Run 24 stress 0.1872714
Run 25 stress 0.3718387
*** Best solution was not repeated -- monoMDS stopping criteria:
20: stress ratio > sratmax
5: scale factor of the gradient < sfgrmin
nmds.2020
Call:
metaMDS(comm = matrix.bray.2020, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2020
Distance: bray
Dimensions: 2
Stress: 0.1522045
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 15 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2020’
stressplot(nmds.2020)
plot(nmds.2020, type="t")
nmds1.2020 <- as.data.frame(scores(nmds.2020, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 4)
)
nmds2.2020 <- as.data.frame(scores(nmds.2020, choices=c(2), display=c("sites")))
nmds_2020_dat <- cbind(nmds1.2020, nmds2.2020) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 4, stroke = 1) +
labs(title = "2020 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 - Grouping by Nutrient Treatment")
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021 <- matrix0 %>%
filter(yr == 2021)
row.2021 <- paste(matrix.2021$nut_trt, matrix.2021$ppt_trt, matrix.2021$yr, matrix.2021$grazing_hist, sep = "_")
matrix.2021 <- matrix.2021[, c(5:79)]
rownames(matrix.2021) <- row.2021
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2021 <- decostand(matrix.2021, "total")
Stress
0.08067789
# relativize data with Bray-Curtis dissimilarity
matrix.2021.1 <- matrix0 %>%
filter(yr == 2021,
grazing_hist == "low")
row.2021.1 <- paste(matrix.2021.1$nut_trt, matrix.2021.1$ppt_trt, matrix.2021.1$yr, matrix.2021.1$grazing_hist, sep = "_")
matrix.2021.1 <- matrix.2021.1[, c(5:79)]
rownames(matrix.2021.1) <- row.2021.1
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2021.1 <- decostand(matrix.2021.1, "total")
nmds.2021.1 <- metaMDS(matrix.bray.2021.1, k=2, trymax = 25)
Run 0 stress 0.08067789
Run 1 stress 0.08067789
... New best solution
... Procrustes: rmse 1.043483e-06 max resid 2.080421e-06
... Similar to previous best
Run 2 stress 0.281194
Run 3 stress 0.08067789
... Procrustes: rmse 7.163003e-06 max resid 1.247917e-05
... Similar to previous best
Run 4 stress 0.08067789
... Procrustes: rmse 1.917534e-06 max resid 2.917769e-06
... Similar to previous best
Run 5 stress 0.08067789
... Procrustes: rmse 1.136837e-06 max resid 1.808009e-06
... Similar to previous best
Run 6 stress 0.1316771
Run 7 stress 0.1393682
Run 8 stress 0.1393365
Run 9 stress 0.1203381
Run 10 stress 0.1407363
Run 11 stress 0.1393679
Run 12 stress 0.203941
Run 13 stress 0.08067789
... New best solution
... Procrustes: rmse 2.590951e-06 max resid 4.031e-06
... Similar to previous best
Run 14 stress 0.3209238
Run 15 stress 0.08067789
... Procrustes: rmse 7.471908e-07 max resid 1.269384e-06
... Similar to previous best
Run 16 stress 0.08067789
... Procrustes: rmse 2.87688e-06 max resid 4.25236e-06
... Similar to previous best
Run 17 stress 0.08067789
... Procrustes: rmse 3.075425e-06 max resid 4.219029e-06
... Similar to previous best
Run 18 stress 0.08067789
... Procrustes: rmse 1.408884e-06 max resid 2.855274e-06
... Similar to previous best
Run 19 stress 0.1203381
Run 20 stress 0.1316771
*** Best solution repeated 5 times
nmds.2021.1
Call:
metaMDS(comm = matrix.bray.2021.1, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2021.1
Distance: bray
Dimensions: 2
Stress: 0.08067789
Stress type 1, weak ties
Best solution was repeated 5 times in 20 tries
The best solution was from try 13 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2021.1’
stressplot(nmds.2021.1)
plot(nmds.2021.1, type="t")
dummy <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites")))
nmds1.2021.1.df <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites"))) %>%
mutate(trt = rownames(dummy)) %>%
separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")
nmds2.2021.1 <- as.data.frame(scores(nmds.2021.1, choices=c(2), display=c("sites")))
nmds_2021_dat.1 <- cbind(nmds1.2021.1.df, nmds2.2021.1) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
# ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
# geom_point(size = 4, stroke = 0.5) +
# labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (LOW) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (LOW) - Grouping by Precipitation Treatment")
Stress
0.08065034
# relativize data with Bray-Curtis dissimilarity
matrix.2021.2 <- matrix0 %>%
filter(yr == 2021,
grazing_hist == "high")
row.2021.2 <- paste(matrix.2021.2$nut_trt, matrix.2021.2$ppt_trt, matrix.2021.2$yr, matrix.2021.2$grazing_hist, sep = "_")
matrix.2021.2 <- matrix.2021.2[, c(5:79)]
rownames(matrix.2021.2) <- row.2021.2
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2021.2 <- decostand(matrix.2021.2, "total")
nmds.2021.2 <- metaMDS(matrix.bray.2021.2, k=2, trymax = 25)
Run 0 stress 0.08065034
Run 1 stress 0.08065038
... Procrustes: rmse 0.0001358001 max resid 0.0002698928
... Similar to previous best
Run 2 stress 0.08065036
... Procrustes: rmse 9.292653e-05 max resid 0.0001874718
... Similar to previous best
Run 3 stress 0.08065038
... Procrustes: rmse 0.0001229219 max resid 0.0002449743
... Similar to previous best
Run 4 stress 0.08065035
... Procrustes: rmse 3.520356e-05 max resid 6.170247e-05
... Similar to previous best
Run 5 stress 0.08065034
... New best solution
... Procrustes: rmse 8.595434e-06 max resid 1.124846e-05
... Similar to previous best
Run 6 stress 0.08554815
Run 7 stress 0.1378309
Run 8 stress 0.08065034
... Procrustes: rmse 5.366069e-05 max resid 9.474028e-05
... Similar to previous best
Run 9 stress 0.08065038
... Procrustes: rmse 0.0001135427 max resid 0.0002111442
... Similar to previous best
Run 10 stress 0.08554818
Run 11 stress 0.08065035
... Procrustes: rmse 3.872519e-05 max resid 7.476398e-05
... Similar to previous best
Run 12 stress 0.08554816
Run 13 stress 0.2444752
Run 14 stress 0.0806506
... Procrustes: rmse 8.06033e-05 max resid 0.0001613536
... Similar to previous best
Run 15 stress 0.08065034
... Procrustes: rmse 2.33919e-05 max resid 4.241008e-05
... Similar to previous best
Run 16 stress 0.08065034
... Procrustes: rmse 3.80112e-05 max resid 8.818892e-05
... Similar to previous best
Run 17 stress 0.08065037
... Procrustes: rmse 0.000114044 max resid 0.0002171751
... Similar to previous best
Run 18 stress 0.1460913
Run 19 stress 0.1378308
Run 20 stress 0.08065034
... Procrustes: rmse 1.981241e-05 max resid 3.849244e-05
... Similar to previous best
*** Best solution repeated 9 times
nmds.2021.2
Call:
metaMDS(comm = matrix.bray.2021.2, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2021.2
Distance: bray
Dimensions: 2
Stress: 0.08065034
Stress type 1, weak ties
Best solution was repeated 9 times in 20 tries
The best solution was from try 5 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2021.2’
stressplot(nmds.2021.2)
plot(nmds.2021.2, type="t")
dummy <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites")))
nmds1.2021.2.df <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites"))) %>%
mutate(trt = rownames(dummy)) %>%
separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")
nmds2.2021.2 <- as.data.frame(scores(nmds.2021.2, choices=c(2), display=c("sites")))
nmds_2021_dat.2 <- cbind(nmds1.2021.2.df, nmds2.2021.2) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (HIGH) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (HIGH) - Grouping by Precipitation Treatment")
Stress: fair
0.1230831
nmds.2021 <- metaMDS(matrix.bray.2021, k=2, trymax = 25)
Run 0 stress 0.1232098
Run 1 stress 0.1688205
Run 2 stress 0.1232098
... Procrustes: rmse 2.53337e-06 max resid 5.886927e-06
... Similar to previous best
Run 3 stress 0.1421227
Run 4 stress 0.1232098
... Procrustes: rmse 4.710548e-06 max resid 9.729523e-06
... Similar to previous best
Run 5 stress 0.1232098
... Procrustes: rmse 1.075313e-05 max resid 2.292137e-05
... Similar to previous best
Run 6 stress 0.1230831
... New best solution
... Procrustes: rmse 0.044648 max resid 0.1321507
Run 7 stress 0.132972
Run 8 stress 0.1507048
Run 9 stress 0.1232098
... Procrustes: rmse 0.04464764 max resid 0.1315101
Run 10 stress 0.1339879
Run 11 stress 0.1230831
... New best solution
... Procrustes: rmse 8.465681e-06 max resid 2.684727e-05
... Similar to previous best
Run 12 stress 0.1232098
... Procrustes: rmse 0.04464354 max resid 0.1315331
Run 13 stress 0.1230831
... New best solution
... Procrustes: rmse 8.230338e-06 max resid 2.599273e-05
... Similar to previous best
Run 14 stress 0.1232098
... Procrustes: rmse 0.04463946 max resid 0.1315179
Run 15 stress 0.132972
Run 16 stress 0.163317
Run 17 stress 0.1230831
... Procrustes: rmse 1.757901e-05 max resid 5.568918e-05
... Similar to previous best
Run 18 stress 0.1339879
Run 19 stress 0.1506988
Run 20 stress 0.1522754
*** Best solution repeated 2 times
nmds.2021
Call:
metaMDS(comm = matrix.bray.2021, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2021
Distance: bray
Dimensions: 2
Stress: 0.1230831
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 13 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2021’
stressplot(nmds.2021)
plot(nmds.2021, type="t")
nmds1.2021 <- as.data.frame(scores(nmds.2021, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 4)
)
nmds2.2021 <- as.data.frame(scores(nmds.2021, choices=c(2), display=c("sites")))
nmds_2021_dat <- cbind(nmds1.2021, nmds2.2021) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 - Grouping by Precipitation Treatment")
Stress: fair
0.1702944
ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2, color = block, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 (Block) - Grouping by Block & Nutrient Treatment")
ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block) - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.low <- matrix0.block %>%
filter(yr == 2021,
block == 3 | block == 4)
row.2021.block.low <- paste(matrix.2021.block.low$nut_trt, matrix.2021.block.low$ppt_trt, matrix.2021.block.low$yr, matrix.2021.block.low$block, sep = "_")
matrix.2021.block.low <- matrix.2021.block.low[, -c(0:4)]
rownames(matrix.2021.block.low) <- row.2021.block.low
Warning: Setting row names on a tibble is deprecated.
matrix.bray.block.2021.low <- decostand(matrix.2021.block.low, "total")
Stress: fair
0.1282296
nmds.2021.block.low <- metaMDS(matrix.bray.block.2021.low, k=2, trymax = 25)
Run 0 stress 0.1282296
Run 1 stress 0.1627532
Run 2 stress 0.1282296
... New best solution
... Procrustes: rmse 2.97563e-05 max resid 6.730899e-05
... Similar to previous best
Run 3 stress 0.1428149
Run 4 stress 0.1282297
... Procrustes: rmse 0.0002046582 max resid 0.0007350626
... Similar to previous best
Run 5 stress 0.214776
Run 6 stress 0.2000888
Run 7 stress 0.1833574
Run 8 stress 0.1282296
... Procrustes: rmse 3.786655e-05 max resid 0.0001314048
... Similar to previous best
Run 9 stress 0.2088749
Run 10 stress 0.1785638
Run 11 stress 0.1785638
Run 12 stress 0.1282296
... Procrustes: rmse 4.723376e-05 max resid 0.0001570042
... Similar to previous best
Run 13 stress 0.2067318
Run 14 stress 0.1428149
Run 15 stress 0.1428149
Run 16 stress 0.1282297
... Procrustes: rmse 0.000151359 max resid 0.0005169016
... Similar to previous best
Run 17 stress 0.1282296
... New best solution
... Procrustes: rmse 2.1128e-05 max resid 6.391618e-05
... Similar to previous best
Run 18 stress 0.1627539
Run 19 stress 0.2006842
Run 20 stress 0.2053314
*** Best solution repeated 1 times
nmds.2021.block.low
Call:
metaMDS(comm = matrix.bray.block.2021.low, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.block.2021.low
Distance: bray
Dimensions: 2
Stress: 0.1282296
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 17 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.block.2021.low’
stressplot(nmds.2021.block.low)
plot(nmds.2021.block.low, type="t")
nmds1.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 3),
block = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 4)
)
nmds2.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(2), display=c("sites")))
nmds_2021_dat.block.low <- cbind(nmds1.2021.block.low, nmds2.2021.block.low) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2, color = block, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 (Block; LOW) - Grouping by Block & Nutrient Treatment")
ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; LOW) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; LOW) - Grouping by Precipitation Treatment")
Stress: fair
0.1412997
nmds.2021.block.high <- metaMDS(matrix.bray.block.2021.high, k=2, trymax = 25)
Run 0 stress 0.1412997
Run 1 stress 0.1856239
Run 2 stress 0.1856239
Run 3 stress 0.1581575
Run 4 stress 0.1581575
Run 5 stress 0.1581575
Run 6 stress 0.1430827
Run 7 stress 0.1412997
... Procrustes: rmse 0.0001557866 max resid 0.000502084
... Similar to previous best
Run 8 stress 0.1412998
... Procrustes: rmse 0.0002316353 max resid 0.0007498162
... Similar to previous best
Run 9 stress 0.1581575
Run 10 stress 0.1581575
Run 11 stress 0.1430823
Run 12 stress 0.1856239
Run 13 stress 0.1581575
Run 14 stress 0.1581576
Run 15 stress 0.1581575
Run 16 stress 0.1856239
Run 17 stress 0.1581575
Run 18 stress 0.1581574
Run 19 stress 0.1430827
Run 20 stress 0.1412998
... Procrustes: rmse 0.0002058184 max resid 0.0006648583
... Similar to previous best
*** Best solution repeated 3 times
nmds.2021.block.high
Call:
metaMDS(comm = matrix.bray.block.2021.high, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.block.2021.high
Distance: bray
Dimensions: 2
Stress: 0.1412997
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.block.2021.high’
stressplot(nmds.2021.block.high)
plot(nmds.2021.block.high, type="t")
nmds1.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 3),
block = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 4)
)
nmds2.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(2), display=c("sites")))
nmds_2021_dat.block.high <- cbind(nmds1.2021.block.high, nmds2.2021.block.high) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, color = block, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 (Block; HIGH) - Grouping by Block & Nutrient Treatment")
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; HIGH) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; HIGH) - Grouping by Precipitation Treatment")
Greenhouse traits
Field traits
“H1: Pairing field results with a greenhouse 15 N experiment will allow us to relate species responses to their ability to uptake organic N”
“In addition, we will conduct a greenhouse experiment to isolate the mechanism (mineralized N vs organic N uptake rates ) by which compost may affect species composition”
Workflow plan:
overlay species trait data as vectors on the NMDS
only greenhouse traits available (only sla data from field)
use PC scores as composite trait axes OR
average the traits for each species
)